warbleR logo

Streamlining acoustic analysis

Grace Smith Vidaurre and Marcelo Araya-Salas
2015-08-30


Here we present a case study of how warbleR functions can be used for bioacoustics, as well as some tips on managing bioacoustic data in R. For more details about function arguments, input or output, please read the documentation for the function in question (e.g. ?querxc).


This vignette can be run without an advanced understanding of R, as long as you know how to run code in you R console. However, knowing more about data and file manipulation would be very helpful in order to modify the code to fit your own needs.


To start, we will use warbleR to download recordings from Xeno-Canto. If you have your own recordings and don’t want to obtain sound files from Xeno-Canto, you can skip the first section and go to Filter recordings by visual inspection.

Prepare for data collection

First, we need to install warbleR

# If this doesn't work, try building warbleR from source package
install.packages("warbleR")
library(warbleR)


And set up a working directory

# Create a new directory
dir.create(file.path(getwd(),"warbleR_example"))
setwd(file.path(getwd(),"warbleR_example"))

# Check the location of the directory
getwd()


Obtain metadata and recordings from Xeno-Canto

Next, we can query the Xeno-Canto database for a species or genus of interest. The function querxc has two types of output:

  1. Metadata of recordings: geographic coordinates, recording quality, recorder, type of signal, etc.

  2. Sound files: Sound files (mp3 format) are returned if the argument download is set to TRUE (default is FALSE).


You can query Xeno-Canto by genus:

# Query Xeno-Canto for all recordings of the hummingbird genus Phaethornis
Phae <- querxc(qword = "Phaethornis", download = FALSE) 
## Obtaining recording list...
## Processing recording information:
## 695 recordings found!
# Find out what kind of metadata we have
names(Phae) 
##  [1] "Recording_ID"      "Genus"             "Specific_epithet" 
##  [4] "Subspecies"        "English_name"      "Recordist"        
##  [7] "Country"           "Locality"          "Latitude"         
## [10] "Longitude"         "Vocalization_type" "Audio_file"       
## [13] "License"           "URL"               "Quality"
View(Phae)


Or you can query by species:

# Query Xeno-Canto for all recordings of the species Phaethornis longirostris
Phae.lon <- querxc(qword = "Phaethornis longirostris", download = FALSE) 
View(Phae.lon)


If you’re interested in the geographic spread of the recording locations, you can use the function xcmaps to visualize locations. xcmaps will create an image file of a map per species in your current directory if img = TRUE. If img = FALSE, maps will be displayed in the RStudio graphic device.

# Image type default is jpeg, but tiff files are often better resolution
# See the documentation for more options if you can't open tiff files
xcmaps(X = Phae, img = TRUE, it = "tiff") 
xcmaps(X = Phae.lon, img = FALSE)

Filter Xeno-Canto recordings by signal type

In most cases, you will need to filter the type of signal you want to analyze. You can filter the recordings prior to the Xeno-Canto download by subsetting the metadata. Then, you can input the filtered metadata back into querxc to download only the selected recordings. There are many ways to filter data in R, and the example below can be modified to fit your own data.

Some of the metadata is not quite consistent across recordings, such as signal type or recording quality. These are characteristics of the recordings that you will need to explore visually with downstream functions before proceeding with data collection and analysis. However, if you are dealing with many recordings, we advise removing the lowest quality recordings (D quality level) or selecting specific vocalization types.

For example, we are interested in assessing the microgeographic variation of long-billed hermit (Phaethornis longirostris) songs. Variation at small geographic scales has been already described in this species (Araya-Salas & Wright 2013). We should proceed by looking for a site with the highest number of songs. Our goal is to search for visible differences in song structure within a site, and then determine whether underlying differences in acoustic parameters are representative of spectrographic distinctiveness.


# Find out number of available recordings
nrow(Phae.lon) 
## [1] 65
# Find out how many types of signal descriptions exist in the Xeno-Canto metadata
levels(Phae.lon$Vocalization_type)
## [1] "300"         "call"        "calls"       "lekking"     "song"       
## [6] "song at lek"
# How many recordings per signal type?
table(Phae.lon$Vocalization_type)
## 
##         300        call       calls     lekking        song song at lek 
##           1           3           1           1          58           1


# There are many levels to the Vocalization_type variable. 
# Some are biologically relevant signals, but most just 
# reflect variation in data entry.

# Luckily, it's very easy to filter the signals we want 
Phae.lon.song <- droplevels(Phae.lon[grep("song", Phae.lon$Vocalization_type, 
                                ignore.case = TRUE), ])

# Check resulting data frame
str(Phae.lon.song) 
## 'data.frame':    59 obs. of  15 variables:
##  $ Recording_ID     : Factor w/ 59 levels "10426","107762",..: 58 57 56 46 37 35 33 30 16 14 ...
##  $ Genus            : Factor w/ 1 level "Phaethornis": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Specific_epithet : Factor w/ 1 level "longirostris": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Subspecies       : Factor w/ 3 levels "","baroni","cephalus": 1 1 2 1 1 1 1 1 1 1 ...
##  $ English_name     : Factor w/ 1 level "Long-billed Hermit": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Recordist        : Factor w/ 7 levels "Fernand DEROUSSEN",..: 6 6 5 3 3 3 3 3 3 3 ...
##  $ Country          : Factor w/ 4 levels "Costa Rica","Ecuador",..: 1 4 2 1 1 1 1 1 1 1 ...
##  $ Locality         : Factor w/ 11 levels "Bonampak, Chiapas",..: 4 9 7 5 5 5 5 5 5 5 ...
##  $ Latitude         : Factor w/ 11 levels "10.3333333333",..: 2 10 7 3 3 3 3 3 3 3 ...
##  $ Longitude        : Factor w/ 10 levels "-78.4651","-79.567",..: 7 4 2 8 8 8 8 8 8 8 ...
##  $ Vocalization_type: Factor w/ 2 levels "song","song at lek": 1 1 2 1 1 1 1 1 1 1 ...
##  $ Audio_file       : Factor w/ 59 levels "http://www.xeno-canto.org/10426/download",..: 58 57 56 46 37 35 33 30 16 14 ...
##  $ License          : Factor w/ 5 levels "http://creativecommons.org/licenses/by-nc-nd/2.5/",..: 3 3 5 4 4 4 4 4 4 4 ...
##  $ URL              : Factor w/ 59 levels "http://www.xeno-canto.org/10426",..: 58 57 56 46 37 35 33 30 16 14 ...
##  $ Quality          : Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1 1 1 1 1 1 ...


# Now, how many recordings per locatity?
table(Phae.lon.song$Locality)
# In case you want more than one signal type:
Phae.lon.sc <- Phae.lon[grep("song|call", Phae.lon$Vocalization_type,ignore.case = TRUE), ]
str(Phae.lon.sc)


# Lets focus on the recordings from La Selva Biological Station
Phae.lon.LS <- Phae.lon.song[grep("La Selva Biological Station, Sarapiqui, Heredia", Phae.lon.song$Locality,
                              ignore.case = FALSE),]

# And only those of the highest quality
Phae.lon.LS <- Phae.lon.LS[Phae.lon.LS$Quality == "A", ]

# We can check if the location coordinates make sense 
# (all recordings should be from a single place in Costa Rica)
# Display a map of these recordings in the graphic device
xcmaps(Phae.lon.LS, img = FALSE)


Once you’re sure you want the recordings, use querxc to download the files. It is also a good idea to save the metadata as .csv files at this point. This data will likely be useful at a later point, especially if you’re aiming for a publication.

# Download sound files
querxc(X = Phae.lon.LS) 

# Save each data frame object as a .csv file 
write.csv(Phae.lon.LS, "Phae_lon.LS.csv", row.names = FALSE)


Convert Xeno-Canto mp3 recordings to wav format

Xeno-Canto maintains recordings in mp3 format, as this compressed format yields smaller files. However, we require wav format for all downstream analyses. Compression from wav to mp3 and back involves information losses, but recordings that have undergone this transformation have been successfully used in research (Medina-Garcia et al. 2015).

To convert mp3 to wav, we can use the warbleR function mp32wav, which relies on a underlying function from tuneR. However, this function does not always work and it remains unclear as to why. This bug should be fixed in future versions of tuneR. If RStudio aborts when running mp32wav, use an mp3 to wav converter online, or download the open source software Audacity (available for Mac, Linux and Windows users). We have made the selected wav files available for download (see next section).

After mp3 files have been converted, we need to check that the wav files are not corrupted and can be read into RStudio (some wav files cannot be read, despite being in wav format).

# Neither of these functions requires arguments
# Always check you're in the right directory beforehand
# getwd()
mp32wav() 

# You can use checkwavs to see if wav files can be read
checkwavs() 

# Let's create a list of all the recordings in the directory
wavs <- list.files(pattern="wav$")

# We will use this list to downsample the wav files so the following analyses go a bit faster
lapply(wavs, function(x) writeWave(downsample(readWave(x), samp.rate = 22050),
                                  filename = x))


### If you were unable to convert mp3 to wav format, you can download the wav files with this code:

# If using Mac, replace "wget" with "curl"
download.file(url = "https://www.dropbox.com/s/wrh5xuidmvn8rno/wavs.RDS?dl=0",
              destfile = "wavs.RDS", method = "wget")

recs <- readRDS(file = "wavs.RDS")

for(i in 1:length(recs))
  writeWave(recs[[i]], filename = names(recs)[i])


Filter recordings by visual inspection

Note: In case you have your own recordings in wav format and have skipped previous sections, you must specify the location of your sound files prior to running downstream functions.

The function lspec is a useful tool for:

  • vocal repertoire analyses
  • filtering by visual inspection

This is the first time we can visualize the recordings since the Xeno-Canto download, and we can make the most of it.

If your research focuses on assessing variation at social or geographic scales, lspec can provide you with important information about how to steer your analysis. If there is an obvious variation in vocalization structure across groups (e.g. treatments or geographic regions), you can focus your analysis on a visual classification of vocalizations. You can use lspec to your advantage here, printing spectrograms on paper and classifying signal types by hand.

Whether or not you decide to proceed with visual classification, lspec allows you to visually inspect the quality of the recording (e.g. amount of background noise) or the type, number, and completeness of the vocalizations of interest. You can discard the image files and recordings that you no longer want to analyze, as this will become very useful for downstream functions.

# Let's first create a subset for playing with arguments 
# This subset is based on the list of wav files we created above
sub <- wavs[c(1,3)]

# ovlp = 10 speeds up process a bit 
# tiff image files are better quality and are faster to produce
lspec(flist = sub, ovlp = 10, it = "tiff")

# We can zoom in on the frequency axis by changing flim, 
# the number of seconds per row, and number of rows
lspec(flist = sub, flim = c(1.5, 11), sxrow = 6, rows = 15, ovlp = 10, it = "tiff")

# Once satisfied with the argument settings we can run all files
lspec(flim = c(1.5, 11), ovlp = 10, sxrow = 6, rows = 15, it = "tiff")

The image files should look like this:

lspec example, make it 700 pxls width

Note that the sound file name and page number are placed in the top right corner.

Now we should inspect the spectrograms. Throwing away image files that are poor quality at first glance (e.g. lots of background noise), will help us in later steps. For instance, the spectrogram for the recording with ID 154143 does not look good, so we will delete those image files. Then, we can select only the sound files that have an image file in the folder. This looks strange for just one recording, but it is really useful when working with many recordings.

# List the image files in the directory
# Change the pattern to "jpeg" if you used that image type
imgs <- list.files(pattern = ".tiff") 

# If the maps we created previously are still there, you can remove them from this list easily
imgs <- imgs[grep("Map", imgs, invert = TRUE)]

# Extract the recording IDs of the files for which image files remain 
kept <- unique(sapply(imgs, function(x){
  strsplit(x, split = "-", fixed = TRUE)[[1]][3]
  }, USE.NAMES = FALSE))

# Now we can get rid of sound files that do not have image files 
snds <- list.files(pattern = ".wav", ignore.case = TRUE) 
file.remove(snds[grep(paste(kept, collapse = "|"), snds, invert = TRUE)])

Data Collection

Automatically detect signals with autodetec

We can move on to data collection with the filtered recordings, using the functions autodetec and manualoc. Both these functions work faster with shorter recordings, but there are ways to deal with larger recordings (an hour long or more).

Here are some points that will help us tailor autodetec for our use:

  1. autodetec has 2 types of output:
    • data frame with recording name, selection, start and end times. These are temporal coordinates that will be passed on to downstream functions when you want to measure acoustic parameters. Save this output as an object, or it will not be saved in the environment, but rather sent to the console.
    • a spectrogram per recording, with red dotted lines marking the start and end of each detected signal, saved in your working directory. In the case of long recordings, as in our example, it’s preferable to create long spectrograms (ls = TRUE). For shorter recordings, spectrograms of the individual selections may work better (ls = FALSE).
  2. Some important detection parameters to fiddle with:
    • threshold controls detection by relative amplitude (%)
    • bp serves as a frequency bandpass filter
    • msmooth controls combination of window length and overlap to smooth signals that have many peaks and would otherwise be detected as multiple signals
    • mindur & maxdur determine the minimum and maximum duration of the signals to be detected

To set these parameters, we need to have some idea of the frequency range and duration of the signals we want to detect. Sectrograms produced above can help us to figure this out. Phaenthornis longirostris songs have frequencies between 2 and 10 kHz and durations between 0.05 and 0.5 s.

If you need to detect all or most of the signals within the recording, play around with different argument values to increase detection accuracy. It may be necessary to do several rounds of optimization with different subsets of your recordings. If just a few signals are needed per recording, a low-accuracy detection could yield enough selections. For instance, if the species you study produces stereotyped signals, like Phaethornis longirostris.

Finally, although autodetec performs automatic signal detection, it doesn’t remove all manual labor from your data collection. Take the time to visually inspect your selections.

# Select a subset of the recordings
wavs <- list.files(pattern = ".wav", ignore.case = TRUE)

# Set a seed so we all have the same results
set.seed(1)
sub <- wavs[sample(1:length(wavs), 3)]

# Run autodetec() on subset of recordings
# Note that no selections are made with these argument settings
autodetec(flist = sub, bp = c(1, 10), threshold = 0, mindur = 0.05, maxdur = 0.5, envt="abs",
          msmooth = c(300, 90), ls = TRUE, res = 100, 
          flim = c(1, 12), wl = 300, set =TRUE, sxrow = 6, rows = 15, 
          redo = FALSE, it = "tiff")

#FIX autodetec TO REDO THE ONES WITH set=T

The image files should look like this (shown below is recording ID 154161):

autodetec image example, make it 700 pxls width

Notice that the images have the argument settings embedded in their file names. This is useful to compare signal detection across combinations of different argument values.

We won’t save the autodetec ouput in an object until we’re satisfied with the detection. To improve our detection we should play around with the arguments values. Below are some detection parameters that work well for these Phaethornis longirotris recordings:

autodetec(flist = sub, bp = c(2, 9), threshold = 20, mindur = 0.09, maxdur = 0.22, 
          envt = "abs", msmooth = c(900,90), ls = TRUE, res = 100, flim = c(1, 12), 
          wl = 300, set =TRUE, sxrow = 6, rows = 15,  redo = TRUE, it = "tiff")

This seems to provide a good detection for most recordings (recording ID 154161):

autodetec image example, make it 700 pxls width

Once we’re satisfied with the detection, we can run the autodetec on all the recordings, removing the argument flist. We also need to save the temporal output at this point.

Phae.ad <- autodetec(bp = c(2, 9), threshold = 20, mindur = 0.09, maxdur = 0.22, 
                     envt = "abs", msmooth = c(900, 90), ls = TRUE, res = 100, 
                     flim= c(1, 12), wl = 300, set =TRUE, sxrow = 6, rows = 15, 
                     redo = TRUE, it = "tiff")

str(Phae.ad)

# Look at the number of selections per sound file 
table(Phae.ad$sound.files)
## 'data.frame':    317 obs. of  5 variables:
##  $ X          : int  1 3 5 7 9 11 13 15 17 19 ...
##  $ sound.files: Factor w/ 6 levels "Phaethornis-longirostris-154070.wav",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ selec      : int  1 3 5 7 9 11 13 15 17 19 ...
##  $ start      : num  0.192 1.466 2.679 3.929 5.624 ...
##  $ end        : num  0.323 1.609 2.838 4.1 5.787 ...
## 
## Phaethornis-longirostris-154070.wav Phaethornis-longirostris-154072.wav 
##                                  17                                 108 
## Phaethornis-longirostris-154123.wav Phaethornis-longirostris-154129.wav 
##                                  30                                  60 
## Phaethornis-longirostris-154138.wav Phaethornis-longirostris-154161.wav 
##                                  22                                  80


Select signals based on signal-to-noise ratio (SNR)


Use snrspecs to prepare signal to noise measurements

Filtering your selected signals by signal to noise ratio (SNR) is often a good idea, but not required. Signals that have a ratio close to 1 (or lower) have very poor quality. A SNR of 1 means the signal and background noise have the same amplitude.

A SNR filter can be applied at any point in your worklow, after using autodetec or manualoc. However, if you just need a sample of the signals in each recording, it would make sense to use the SNR functions to perform another quality filter prior to making acoustic measurements. Like the other functions downstream of autodetec or manualoc, the signal to noise functions require the start and end time of the signals.

snrspecs is another function in the family of spectrogram creators. It has very similar arguments to specreator, but it also has additional arguments for picking a margin over which to measure noise. These margins are very important for calculating SNR, especially when you’re measuring signals with short silence in between. You want to be sure to pick a noise margin that doesn’t overlap neighboring signals.

# A margin that's too large causes other signals to be included in the noise measurement
# Re-initialize X as needed, for either autodetec or manualoc output

# Let's try it on 10% of the selections so it goes a faster
# Set a seed first, so we all have the same results
set.seed(5)

X <- Phae.ad[sample(1:nrow(Phae.ad),(nrow(Phae.ad)*0.1)), ]

snrspecs(X = X, flim = c(2, 11), snrmar = 0.5, mar = 0.7, it = "tiff")

The image files should look like this:

snrpecs

This margin is far too large! It’s overlapping the whole signal. Let’s try with smaller margins.

# This smaller margin is better
snrspecs(X = X, flim = c(2, 11), snrmar = 0.2, mar = 0.7, it = "tiff")

snrpecs

These smaller margins seem to work well. Now we can run the function over all the selections and inspect all spectrograms to be sure that the margin(s) still hold.

snrspecs(X = Phae.ad, flim = c(2, 11), snrmar = 0.2, mar = 0.7, it = "tiff")


Calculate signal to noise ratio for recordings

Once you’ve picked a margin for all recordings, you can move forward with the SNR calculation. This calculation can allow you to remove poor quality recordings with a SNR close to 1. Since you’ve already performed several visual filters in the workflow, this step often isn’t necessary, but it can provide you with quantitative information about recording quality.

We will measure SNR on every other selection to speed up the process

Phae.snr <- sig2noise(X = Phae.ad[seq(1, nrow(Phae.ad), 2), ], mar = 0.04)

As we just need a few songs to characterize each sound file and individual, we can select selections with the highest SNR per sound file. In this example, we will choose 5 selections with the highest SNRs.

Phae.hisnr <- Phae.snr[ave(-Phae.snr$SNR, Phae.snr$sound.files, FUN = rank) <= 5, ]

# Double check the number of selection per sound files 
table(Phae.hisnr$sound.files)
## 
## Phaethornis-longirostris-154070.wav Phaethornis-longirostris-154072.wav 
##                                   5                                   5 
## Phaethornis-longirostris-154123.wav Phaethornis-longirostris-154129.wav 
##                                   5                                   5 
## Phaethornis-longirostris-154138.wav Phaethornis-longirostris-154161.wav 
##                                   5                                   5

At this point would be a good idea to save the selections as a file

write.csv(Phae.hisnr, "Phae_lon_autodetec_selecs.csv", row.names = FALSE)


Manually select signals with manualoc

In some cases manual selection may be preferable, especially if you have shorter recordings or if the automatic detection is not as accurate as you’d like.

manualoc is a function that provides a graphical interface in which you can select the start and end of the signals. It can often run slowly, depending on the size of the sound files. However, manualoc can be very useful when you can get away with selecting only a few signals, perhaps if your species has stereotyped signals, or if you’re interested in a specific element of a complex signal across recordings.

We recommend reading the documentation for manualoc prior to running this example. Once you’ve done so, here are some points to keep in mind:

  1. The sole output for this function is a .csv file:
    • contains the time coordinates, selection information and any comments made
    • similar to the autodetec output, these coordinates will be used in downstream functions
  2. Be very precise with your clicks
    • stray clicks will cause manualoc to fail
    • don’t double-click, instead click once and wait for blue bubble to pop up before clicking again
    • suspect a stray click when manualoc starts responding to single clicks
    • if so, stop manualoc with Stop button, or with red Stop button in RStudio console
    • manualoc retains all previous selections in the .csv file and will start up where you left off
  3. Selections can be deleted directly in the manualoc interface
    • use Del-sel button
  4. Selections can also be deleted in the manualoc_output.csv
    • stop manualoc, open the .csv
    • delete the rows corresponding to the unwanted selection(s) prior to starting manualoc again
  5. Run manualoc within the expected frequency range for your species
    • use argument flim to facilitate signal selection
  6. Run manualoc with oscillograms enabled to improve signal selection
    • when osci = TRUE, the oscillogram or waveform serve as a visual aid
    • use changes in amplitude to select start and end of the signal
    • oscillogram will print to screen when the resolution of the projected spectrogram
      improves (depends on the seltime argument)
    • if seltime = 2, the oscillogram will show up for selections <= 2 seconds

Some other uses for manualoc:

  1. manualoc can be used in combination with autodetec if you have large recordings:
    • each file will take a long time to load
    • but you can select specific time points to break up the recording
    • then you can feed these time coordinates to autodetec using the data frame argument X
    • this speeds up the automated detection process
    • this can help customize autodetec if you have recordings with different noise or playback treatments
  2. manualoc can also be used for visual classification
  • if you don’t have a printer to print and mark long spectrograms by hand
  • run manualoc with selcomm = TRUE
  • mark individual selections with song or element types using selcomm
  • use specreator to create spectrograms with selcomm text and check visual classifications

Note that you can stop the function at any point by clicking twice on the stop button.

# Run manualoc() with frequency range set for Phaethornis longirostris
# Recording comments are enabled to mark recording quality
# Selection comments enabled to include visual classifications
manualoc(flim = c(1, 11), reccomm = TRUE, selcomm = TRUE, osci = TRUE)

# Read manualoc() output back into RStudio as an object
# This data frame object can be used as input for the functions that follow
manualoc_out <- read.csv("manualoc_output.csv", header = TRUE)

The graphic device should display something like this

manualoc view


Visualize autodetec or manualoc selections with specreator

specreator serves as yet another option for visual inspection, although at the level of individual selections made through autodetec and manualoc.

Like the other members of the spectrogram-creating family, specreator contains many options related to graphical parameters. With some fiddling around, it’s possible to make images of publication quality. However, some of these graphical parameters do not play well together (especially osci, gr, sc), see the documentation for suggestions.

# Create a subset of 5 recordings analyzed by autodetec() or manualoc()
# Speeds up process of playing around with arguments 
# Run either line below to reinitialize X with either autodetec 
# or manualoc subset as desired

set.seed(50)
X <- Phae.hisnr[sample(1:nrow(Phae.hisnr), 5), ]

# Plot selection lines from manualoc() or autodetec()
specreator(X, osci = FALSE, line = TRUE, wl = 300, flim = c(1, 11), it = "tiff")

# Change frequency limits of y-axis
specreator(X, flim = c(1, 11), osci = TRUE, line = TRUE, wl = 300, it = "tiff")

# Change width of spectrogram to be proportional to signal duration
specreator(X, flim = c(1, 11), osci = TRUE, line = TRUE, propwidth = TRUE, wl = 300, it = "tiff")

# Change spectrogram size 
# Changing inner.mar and outer.mar arguments improves picsize results
specreator(X, flim = c(1, 11), osci = TRUE, line = TRUE, picsize = 1.5, wl = 300, 
           ovlp = 90, inner.mar = c(4, 4.5, 2, 1), outer.mar = c(4, 2, 2, 1), it = "tiff")

# Run function for all recordings, with final argument settings
specreator(Phae.hisnr, flim = c(1, 11), osci = TRUE, line = TRUE, wl = 300,
           ovlp = 90, it = "tiff", res = 300)

specreator example

Measure acoustic parameters

Visualize frequency measurements with trackfreqs

Prior to calculating acoustic measurements, it’s good practice to visualize the accuracy of some important measurements, namely frequency measurements. The function trackfreqs is the last in the family of spectrogram-creators. It allows you to create spectrograms with dominant frequency and fundamental frequency measurements plotted on top of each selected signal.

In general, the fundamental frequency measurements are not as reliable as the dominant frequency measurements. When aocustic measurements are performed in specan, the fundamental frequency reported is the mean of all the individual fundamental frequency measurements, and is usually more accurate. Use trackfreqs on all the recordings for which you want to measure acoustic parameters. Scroll through all the spectrograms to get a feeling for how well the frequency measurements will be performed across your recordings.

Like it’s sister functions, trackfreqs has many graphical arguments. It has additional graphical arguments to change colors of the plotting symbols, and size and position of legend labels. These arguments will largely depend on the nature of your selections.

# Note that the dominant frequency measurements are almost always more accurate
trackfreqs(Phae.hisnr, flim = c(1, 11), bp = c(2, 12))

# Play around with the colors and sizes of the symbols
# see par() and points() in RStudio help for more details
trackfreqs(Phae.hisnr, flim = c(1, 11), bp = c(6, 8), col = c("purple", "orange"),
           pch = c(17, 3), res = 300)

# We can change the lower end of bandpass to make the frequency measurements more precise
# If the frequency measurements look acceptable with this bandpass setting,
# that's the setting we should use when running specan() 
trackfreqs(Phae.hisnr, flim = c(1, 11), bp = c(2, 12), col = c("purple", "orange"),
           pch = c(17, 3), res = 300)

trackfreqs image

As fundamental frequency does not seem to be adequately tracked, we will remove it from the acoustic parameters measured in the next step.


Batch-process acoustic measurements with specan

We’re close to finshing the warbleR workflow. We can now perform acoustic measurements with the function specan. This function calculates 22 acoustic parameters across all the specified recordings. It’s a batch process that is much faster than calculating measurements one recording at a time. specan uses and customizes several functions available in the package seewave.

If you require more customized acoustic measurements, you can add your own customizations of seewave functions to the specan code. We suggest that you create your own copy of specan, rather than modifying the package version. Changing the specan code will require getting used to trouble-shooting and debugging in R, unless you’re already comfortable writing your own functions.

# specan() uses the time coordinates in the autodetec or manualoc output
# It will measure acoustic parameters within the start and end times per selection

# Use the bandpass filter to your advantage, to filter out low or high background
# noise before performing measurements
# The amplitude threshold will change the amplitude at which noises are
# detected for measurements 
params <- specan(Phae.hisnr, bp = c(1, 11), threshold = 15)

View(params)

str(params)

# As always, it's a good idea to write .csv files to your working directory
## 'data.frame':    30 obs. of  24 variables:
##  $ sound.files: Factor w/ 6 levels "Phaethornis-longirostris-154070.wav",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ selec      : int  13 15 19 27 29 104 126 132 144 148 ...
##  $ duration   : num  0.176 0.176 0.176 0.159 0.163 ...
##  $ meanfreq   : num  6.53 6.5 6.49 6.61 6.51 ...
##  $ sd         : num  1.77 1.78 1.74 1.84 1.78 ...
##  $ median     : num  6.45 6.34 6.34 6.58 6.42 ...
##  $ Q25        : num  5.51 5.5 5.46 5.39 5.44 ...
##  $ Q75        : num  7.66 7.66 7.61 7.78 7.61 ...
##  $ IQR        : num  2.15 2.16 2.15 2.39 2.17 ...
##  $ skew       : num  1.96 2.26 2.21 2.18 2.04 ...
##  $ kurt       : num  7.73 9.03 9.28 10.05 8.71 ...
##  $ sp.ent     : num  0.943 0.941 0.94 0.947 0.944 ...
##  $ sfm        : num  0.594 0.613 0.572 0.626 0.603 ...
##  $ mode       : num  5.87 6.12 5.88 6.11 6.13 ...
##  $ centroid   : num  6.53 6.5 6.49 6.61 6.51 ...
##  $ peakf      : num  5.84 5.6868 5.7776 6.6079 0.0427 ...
##  $ meanfun    : num  4.04 4.38 4.91 4.81 4.83 ...
##  $ minfun     : num  0.286 0.882 2.756 2.756 2.756 ...
##  $ maxfun     : num  7.35 7.35 7.35 7.35 7.35 7.35 7.35 7.35 7.35 7.35 ...
##  $ meandom    : num  5.9 6.52 6.2 6.87 6.79 ...
##  $ mindom     : num  5.17 5.86 5.34 6.12 5.9 ...
##  $ maxdom     : num  7.19 7.75 7.71 8.01 8.05 ...
##  $ dfrange    : num  2.02 1.89 2.37 1.89 2.15 ...
##  $ modindx    : num  0.447 0.576 0.348 0.627 0.457 ...

Now let’s remove parameters derived from fundamental frequency (based on trackfreqs results)

params <- params[, grep("fun|peakf", colnames(params), invert = TRUE)]


Analysis of geographic variation using specan measurements

Now we can evaluate whether the observed variation in song structure is actually reflected by the acoustic parameters we just measured. For this we will conduct a Principal Component Analysis on scaled (z-transformed) acoustic parameters and look at the grouping of songs (data points) in the scatter plot.

# Run the PCA with only numeric variables of params
pca <- prcomp(x = params[, sapply(params, is.numeric)], scale. = TRUE)

# Check loadings
summary(pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4    PC5    PC6     PC7
## Standard deviation     2.9413 2.2746 1.35134 1.08980 0.9389 0.6195 0.56543
## Proportion of Variance 0.4553 0.2723 0.09611 0.06251 0.0464 0.0202 0.01683
## Cumulative Proportion  0.4553 0.7276 0.82375 0.88626 0.9327 0.9528 0.96968
##                            PC8     PC9    PC10    PC11   PC12    PC13
## Standard deviation     0.52419 0.32400 0.28518 0.24237 0.1744 0.11011
## Proportion of Variance 0.01446 0.00553 0.00428 0.00309 0.0016 0.00064
## Cumulative Proportion  0.98414 0.98967 0.99395 0.99704 0.9986 0.99928
##                           PC14    PC15    PC16      PC17      PC18
## Standard deviation     0.08364 0.06496 0.04952 2.748e-15 3.341e-16
## Proportion of Variance 0.00037 0.00022 0.00013 0.000e+00 0.000e+00
## Cumulative Proportion  0.99965 0.99987 1.00000 1.000e+00 1.000e+00
##                             PC19
## Standard deviation     1.392e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00
# Extract PCA scores
pcascor <- as.data.frame(pca[[5]])

# Plot the 2 first PCs
plot(pcascor[, 1], pcascor[, 2], col = as.numeric(params$sound.files), pch = 20, 
     cex = 3, xlab = "PC1", ylab = "PC2")

# Add recordings/individuals labels 
x <- tapply(pcascor[, 1], params$sound.files, mean)
y <- tapply(pcascor[, 2], params$sound.files, mean)

labs <- gsub(".wav", "", unique(sapply(as.character(params$sound.files), function(x){
  strsplit(x, split = "-", fixed = TRUE)[[1]][3]
  }, USE.NAMES = FALSE)))

text(x, y, labs)

It seems like the songs are grouped by sound file, or individual signature. Let’s look at the song type level. First, we need to classify the songs by song type. We can check the spectrograms we previously created to do this.

song types, make it 700 pxl width

Songs from sound files 154070 and 154072 seem to belong to the same song type. Sound files 154129 and 154161 represent a different song type. Finally, the songs from each of the other 2 sound files has a unique structure, so each one represents a different song type. We can add this information to the plot by using symbols to represent song types.

# Create a song type variable

# First, extract recording ID
songtype <- gsub(".wav", "", sapply(as.character(params$sound.files), function(x){
  strsplit(x, split = "-", fixed = TRUE)[[1]][3]
  }, USE.NAMES = FALSE))

# Now change IDs for letters representing song types
songtype <- gsub("154070|154072", "A", songtype)
songtype <- gsub("154129|154161", "B", songtype)
songtype <- gsub("154123", "C", songtype)
songtype <- gsub("154138", "D", songtype)

# Add song type as a variable representing symbol type
plot(pcascor[, 1], pcascor[, 2], col = as.numeric(params$sound.files), pch = as.numeric(as.factor(songtype)), 
     cex = 3, xlab = "PC1", ylab = "PC2")

# Add song type labels 
x <- tapply(pcascor[, 1], songtype, mean)
y <- tapply(pcascor[, 2], songtype, mean)

text(x, y, unique(songtype), cex = 1.5)

Songs of the same song type are more similar (they cluster together). This PCA confirms that the visually distinctive song types are representative of underlying acoustic difference. Likewise, it also confirms that we can detect biologically relevant differences from sound files that have undergone mp3 compression and conversion back to wav format.


References

Araya-Salas, M. and T. Wright. 2013. Open-ended song learning in a hummingbird. Biology Letters. 9 (5). doi: 10.1098/rsbl.2013.0625 PDF

Medina‐García, Angela, M. Araya‐Salas, and T. Wright. 2015. Does vocal learning accelerate acoustic diversification? Evolution of contact calls in Neotropical parrots. Journal of Evolutionary Biology. doi: 10.1111/jeb.12694 PDF